;;........................................................................ 
;; Easy Aerosol diag
;;
;; Clear-sky shortwave flux at TOA
;; Clear-sky shortwave flux at surface 
;; Clear-sky shortwave flux in atmosphere
;;........................................................................ 
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"   

begin

  hres = "T42" 
  yy = "1990-1999" 
  yy = "1980-1999" 
  yy = "1980-1989" 
  yy = "1979-2008" 
  tp = "amip" 

;;  hres = getenv("hres") 
;;  yy = getenv("yy") 
;;  tp = getenv("tp") 

  xx = "_" 

  plotname = "./figure/figure_clearsky_zonal_vs_3plume_line" + xx + yy + xx + hres + xx + tp 
  plotname = "figure_clearsky_zonal_vs_3plume_line" + xx + yy + xx + hres + xx + tp 

  myfont = 0.025

  if(tp.eq."amip") then 
     cs1 = "EXP04"  
     cs2 = "EXP05"  
     cs3 = "EXP06"  
  else 
     cs1 = "EXP01"  
     cs2 = "EXP02"  
     cs3 = "EXP03"  
  end if 

  fna = cs1+"_"+hres+xx+yy+".nc"
  fnb = cs2+"_"+hres+xx+yy+".nc" 
  fnc = cs3+"_"+hres+xx+yy+".nc" 

  print("read file : "+fna) 
  print("read file : "+fnb) 
  print("read file : "+fnc) 

  system("ls -l "+fna) 
  system("ls -l "+fnb) 
  system("ls -l "+fnc) 

  fla = addfile(fna,"r") 
  flb = addfile(fnb,"r") 
  flc = addfile(fnc,"r") 

  lat  = fla->lat
  lon  = fla->lon
  time = fla->time

  ;; TOA 

  atoa = fla->FSNTC
  btoa = flb->FSNTC
  ctoa = flc->FSNTC

  pzn_toa = atoa
  pzn_toa = btoa - atoa

  p3p_toa = atoa
  p3p_toa = ctoa - atoa

  zonal_pzn_toa = dim_avg(pzn_toa) 
  zonal_p3p_toa = dim_avg(p3p_toa) 

  ;; Surface 

  asfc = fla->FSNSC
  bsfc = flb->FSNSC
  csfc = flc->FSNSC

  pzn_sfc = asfc
  pzn_sfc = bsfc - asfc

  p3p_sfc = asfc
  p3p_sfc = csfc - asfc

  zonal_pzn_sfc = dim_avg(pzn_sfc) 
  zonal_p3p_sfc = dim_avg(p3p_sfc) 

  ;; Atmosphere

  aatm = atoa 
  batm = atoa 
  catm = atoa 

  aatm = atoa - asfc 
  batm = btoa - bsfc 
  catm = ctoa - csfc 

  pzn_atm = aatm
  pzn_atm = batm - aatm

  p3p_atm = aatm
  p3p_atm = catm - aatm

  zonal_pzn_atm = dim_avg(pzn_atm) 
  zonal_p3p_atm = dim_avg(p3p_atm) 

  ;; Surface

  aaod = fla->PSA_TAU 
  baod = flb->PSA_TAU 
  caod = flc->PSA_TAU 

  pzn_aod = aaod
  pzn_aod = baod - aaod

  p3p_aod = aaod
  p3p_aod = caod - aaod

  zonal_pzn_aod = dim_avg(pzn_aod)
  zonal_p3p_aod = dim_avg(p3p_aod)

;; test 
;; pzn_atm = pzn_atm / pzn_aod * p3p_aod 


  nlon = dimsizes(lon) 
  nlat = dimsizes(lat) 

  ;;................................
  ;; plotting 
  ;;................................

  plot    = new ( 4, "graphic")

  wks  = gsn_open_wks ("pdf",plotname) 

  colormap = "ViBlGrWhYeOrRe"
  colormap = "testcmap"
  colormap = "amwg"
  colormap = "BlueYellowRed" 

  gsn_define_colormap(wks,colormap)         ; choose a colormap


  res  = True
  res@gsnDraw                = False
  res@gsnFrame               = False
  ;;res@tiYAxisString          = "PDF (%)"
  ;;res@tiXAxisString          = "10m wind speed (m s~S~-1~N~)"

  res@tiMainOn         = False   ; no title
  res@cnFillOn         = True    ; turn on color
  res@cnLinesOn        = False   ; turn off contour lines
  res@cnLineThicknessF = 1.     ; thicker lines
  ;;res@cnLevelSelectionMode  = "ManualLevels"  ; set manual cn levels
  ;;res@cnMinLevelValF   =    0.     ; set min contour level
  ;;res@cnMaxLevelValF   = 2.     ; set max contour level
  ;;res@cnLevelSpacingF  = 0.1
  res@cnInfoLabelOn    = False    ; no contour labels
;;  res@cnFillMode       = "RasterFill"  ; turn on raster mode
 ;res@cnFillDrawOrder  = "PreDraw"     ; draw contours before continents

  res@gsnSpreadColors = True    ; use full colormap



;;----------------------------------------------------------------------
;; lable bar
;;----------------------------------------------------------------------

  res@lbLabelBarOn       = True       ; no individual label bar
  res@lbLabelStride      = 1          ; every other label bar label
  res@lbOrientation      = "Vertical" ;;Horizontal" ; vertical label bar
  res@lbLabelFontHeightF = myfont * 0.8
  res@pmLabelBarWidthF   = 0.08
  res@pmLabelBarHeightF  = 0.4
  ;;res@pmLabelBarOrthogonalPosF =  0.18


  ;;colors = (/"white","black","White","RoyalBlue","LightSkyBlue",\
  ;;            "PowderBlue","lightseagreen","PaleGreen","Wheat","Brown",\
  ;;            "Pink","darkgreen","darkorange","red","grey"/)

  ;;;;res@xyLineThicknessF       = 2
  ;;res@xyLineThicknesses = 4. + 0.*fspan(0,10,10)
  ;;;;res@xyDashPatterns    = (/0.,1.,2.,0.,1.,2./) ;;0.*fspan(0,10,10) 
  ;;res@xyDashPatterns    = (/0.,1.,0.,1./) ;;0.*fspan(0,10,10) 
  ;;;;res@xyLineColors      = (/"RoyalBlue","grey","Wheat","blue","red","darkgreen","darkorange","brown","lightseagreen","black"/)    ; change line color
  ;;;;res@xyLineColors      = (/"blue","blue","green","red","red","orange"/) 
  ;;res@xyLineColors      = (/"black","red","red","red"/) 

  res@gsnStringFontHeightF     = myfont
  res@tmXBLabelFontHeightF = myfont * 0.8 
  res@tmYLLabelFontHeightF = myfont* 0.8 

  res@tiXAxisFontHeightF = myfont
  res@tiYAxisFontHeightF = myfont

  res@cnLevelSelectionMode  = "ExplicitLevels"  ; set manual cn levels
  res@cnLevels = (/-8,-6,-4,-2,-1,0,1,2,4,6,8/) 

  res@trYMinF  = -6.0                   ; min value on y-axis
  res@trYMaxF  =  3.7                   ; max value on y-axis


  res@xyLineThicknesses = (/3.0,3.0/)               ; make 2nd lines thicker


  res@tiYAxisString    = "~F33~D~F~ TOA clear-sky SW "  
  ;;res@gsnLeftString = "a) ~F33~D~F~ TOA clear-sky SW flux perturbation" 
  res@gsnRightString = "W m~S~-2~"
  res@gsnCenterString        = "" 


  res@vpHeightF= 0.4                    ; change aspect ratio of plot
  res@vpWidthF = 0.8                  

  res@trYMinF  = -5.0                   ; min value on y-axis
  res@trYMaxF  =  3.0                   ; max value on y-axis


  res@xyLineColors      = (/"blue"/)          ; change line color
  plot(0) = gsn_csm_xy(wks, lat, zonal_pzn_toa(0,:), res)
  res@xyLineColors      = (/"red"/)          ; change line color
  plot0 = gsn_csm_xy(wks, lat, zonal_p3p_toa(0,:), res)
  overlay(plot(0),plot0) 

  res@trYMinF  = -8.0                   ; min value on y-axis
  res@trYMaxF  =  2.0                   ; max value on y-axis

  ;;res@gsnLeftString = "b) Surface clear-sky SW flux perturbation" 
  res@tiYAxisString    = "~F33~D~F~ SFC clear-sky SW " 
  res@gsnRightString = "W m~S~-2~"
  res@gsnCenterString        = "" 
  res@xyLineColors      = (/"blue"/)          ; change line color
  plot(1) = gsn_csm_xy(wks, lat, zonal_pzn_sfc(0,:), res)
  res@xyLineColors      = (/"red"/)          ; change line color
  plot1 = gsn_csm_xy(wks, lat, zonal_p3p_sfc(0,:), res)
  overlay(plot(1),plot1) 

  res@trYMinF  = -1.0                   ; min value on y-axis
  res@trYMaxF  =  5.0                   ; max value on y-axis

  ;;delete(res@cnLevels) 
  ;;res@gsnLeftString = "c) ATM clear-sky SW flux perturbation" 
  res@tiYAxisString    = "~F33~D~F~ ATM clear-sky SW "
  res@gsnRightString = "W m~S~-2~"
  res@gsnCenterString        = "" 
  res@xyLineColors      = (/"blue"/)          ; change line color
  plot(2) = gsn_csm_xy(wks, lat, zonal_pzn_atm(0,:), res)
  res@xyLineColors      = (/"red"/)          ; change line color
  plot2 = gsn_csm_xy(wks, lat, zonal_p3p_atm(0,:), res)
  overlay(plot(2),plot2) 

  res@trYMinF  =  0.0                   ; min value on y-axis
  res@trYMaxF  =  0.1

  res@tiYAxisString    = "AOD"
  res@gsnRightString = ""
  res@gsnCenterString        = ""
  res@xyLineColors      = (/"blue"/)          ; change line color
  plot(3) = gsn_csm_xy(wks, lat, zonal_pzn_aod(0,:), res)
  res@xyLineColors      = (/"red"/)          ; change line color
  plot3 = gsn_csm_xy(wks, lat, zonal_p3p_aod(0,:), res)
  overlay(plot(3),plot3)


  resP    = True
  resP@txString  = ""
  resP@gsnPanelYWhiteSpacePercent = 5
  resP@gsnPanelXWhiteSpacePercent = 5
  gsn_panel(wks,plot,(/2,2/),resP)

  print("plotname : "+plotname) 

end







